;---------------------------------------------------------------------
; calculate radon concentration at transcom sites 
; Kai Zhang, MPI-M 
;---------------------------------------------------------------------
  load "$KAILIB_ROOT/calllib.ncl" 
  load "$KAILIB_ROOT/callprt.ncl" 
  load "$KAILIB_ROOT/callvar.ncl" 
  load "$KAILIB_ROOT/lib_interp.ncl" 
  load "$KAILIB_ROOT/lib_radon.ncl" 
  load "$KAILIB_ROOT/lib_transcom.ncl" 

begin 

  cmd = printcls(5) 
  cmd = printmpi 
  cmd = printcls(5) 

  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  ;; load.ncl gives following information 
  ;; generated by csh shell  
  ;;   xlon 
  ;;   xlat
  ;;   xhgt
  ;;   filenama 
  ;;   filenamo 
  ;;   scale 
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  load "load.ncl" 

  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  ;; get site information from transcom data file
  ;; data read from lib_transcom  
  ;;   lat
  ;;   lon
  ;;   height 
  ;;   y/n land site
  ;;   y/n above sea level 
  ;;   site id 
  ;;   number of site  
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  alat = transcom_lat(0)
  alon = transcom_lon(0)
  ahgt = transcom_hgt(0)
  alnd = transcom_land(0)
  aslv = transcom_absl(0)
  anum = transcom_id(0)
  ns   = transcom_ns(0)

  cmd = printcls(3)
  print(sprinti("%0.3i",anum)+"   "+ \
        sprintf("%8.3f",alat)+"   "+ \
        sprintf("%8.3f",alon)+"   "+ \
        sprintf("%8.3f",ahgt)+"   "+ \
        sprinti("%0.1i",alnd)+"   "+ \
        sprinti("%0.1i",aslv)        \
       )
  cmd = printcls(3)

  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  ;; read data  
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  fila = addfile(filenama,"r") 

  print("") 
  print("input  filename: "+filenama) 
  print("output filename: "+filenamo) 
  print("") 

  radon = fila->RN222
  lon   = fila->lon
  lat   = fila->lat
  time  = fila->time

  nt = dimsizes(time) 

  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  ;; define output data format 
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  rna = new((/nt,ns/),"float") 
  rna!0 = "time" 
  rna&time = time
  rna!1 = "site" 
  rna&site = anum 

  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  ;; interpolation 
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  ;rnt = (/ linint2_points (lon,lat(::-1),radon(:,:,::-1,:),True,alon,alat,0) /) 
  ; 
  ;do is = 0, ns-1 
  ;   xlon = alon(is) 
  ;   xlat = alat(is) 
  ;   xhgt = ahgt(is) 
  ;   ilev = location_level(xlon,xlat,xhgt,1,0)  
  ;   rna(:,is) = (/ rnt(:,ilev,is) /) 
  ;end do 

  do is = 0, ns-1 
     xlon = alon(is) 
     xlat = alat(is) 
     xhgt = ahgt(is) 
     ilev = location_level(xlon,xlat,xhgt,1,0)  
     rna(:,is) = (/ linint2_points(lon,lat(::-1),radon(:,ilev,::-1,:),True,xlon,xlat,0) /) 
  end do 

  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  ;; unit conversion using lib_radon  
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  rna = RN222_mBqm3(rna) 

  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  ;; scale it if needed 
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  rna = rna * scale 

  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  ;; check data 
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  cmd = printcls(3)
  print(sprinti("%0.3i",anum)+"   "+ \
        sprintf("%8.3f",alat)+"   "+ \
        sprintf("%8.3f",alon)+"   "+ \
        sprintf("%8.3f",ahgt)+"   "+ \
        sprinti("%0.1i",alnd)+"   "+ \
        sprinti("%0.1i",aslv)+"   "+ \
        sprintf("%9.3f",xdimavg(rna,"time"))  \
       )
  cmd = printcls(3)

  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  ;; output data 
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  system("rm "+filenamo) 
  filo = addfile(filenamo,"c") 
  filo->radon=rna

  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  ;; program end 
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  cmd = printbye

end 

